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The Isinig model is a famous and applicable model in the statistical 
mechanics. In this project, I propose to consider the Ising model and its 
computational simulation. In this text, I try to review the main idea of 
the Ising model and Metropolis algorithm, the computational approach to 
this problem. Finally, I report the result of my own simulation and 
compare its results with analytical solution as a conclusion. 
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INTRODUCTION 

Statistical mechanics is a very active area of research and there are many open problems. In 
fact, the general phenomenon of phase transitions is still actively being researched. An 
example of these phenomena is ferromagnetic transition in the Curie temperature. In this 
phenomenon, Ferromagnetic materials, when heated, eventually lose their magnetic 
properties. This loss becomes complete above the Curie temperature, named after the French 
physicist Pierre Curie, who discovered it in 1895 [1-2]. 

Paramagnetism is the normal induced magnetization of a material when it is put into a 
magnetic field. Paramagnetism could be explained by assuming that all spins are independent. 
Ferromagnetism is the spontaneous magnetization occurring in certain metals such as iron and 
nickel. To explain this phenomenon we must take into account the interaction between spins 
[2]. 

The Ising model is an attempt to simulate the structure of this phase transition. Originally, it 
was invented for the phase transition of ferromagnets at the Curie temperature; however, in 
the course of time it was realized that with only slight changes the model can also be applied 
to other phase transitions, like order-disorder transitions in binary alloys. Furthermore, the 
model may be applied to several modem problems of many particle physics, for instance for 
the description of so-called spin glasses. These are metals having amorphous instead of 
crystalline structures, which have the interesting property of no vanishing entropy at T = 0. 
Recently, it has been realized that Ising's idea (in modified form) could also explain pattern 
recognition in schematic neural networks. Thus, this model gains more and more importance 
for the development of models for the human brain [3-8]. 

The main idea of Ising model is one orientation spin-interaction of fixed point on the lattice. 
This simplification is based on quantum mechanics concepts 'exchange force' [1]. It's 
meaning that the only outer electrons (conduction electrons) on the matter have interaction 
and their interaction is only in z-direction. On the other hand, this means that the characters of 
lattice like dimension, size, morphology and structure are variables of this model. We have to 
resume these quantities to have a specific sample of Ising model. 

The other consideration of Ising model is that the only nearest neighbors interactions are 
necessary and others can be neglect [5]. The finally assumption is that the lattice in statistical 
limit are periodic [6]. 
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Figure 1. The 2-demension lattice; 
(a) Square (4 nearest neighbor); (b) Triangular (6 nearest neighbor) (c) Hexagunal (3 nearest neighbor). [7] 
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There are exact analytical solutions for 1 and 2 dimensions Ising models. In this project I have 
to simulate the 2-dimension Ising model with Square lattice. To see exact solution of this 
structure with addition theoretical detail, refer to Riechl [7]. 



Figure 2. The 2-demension Ising with square lattice. 
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Figure 3. The summary of analytical and experimental results T c = 2.269 — . 



COMPUTATIONAL APPROACH 

The computational simulation of Ising model is based on the two facts: The statistical nature 
of events and the Ising Hamiltonian [9]. From statistical mechanics, we know the probability 
of every event A is proportional to Boltzmann factor: 

P04)ocexp(-^;) 
so the transition probability of A -> B phase transition is: 

As 

P(7l^£)ocexp(-— ) 

The computational algorithm which applies this factor at a problem is Metropolis algorithm 
[9], which I needed to simulate the Ising model. In this section, I will introduce the 
Metropolis algorithm and then present my own simulation algorithm. 



(a) Metropolis Algorithm Steps [9]: 

1. Start with an arbitrary spin configuration (initial conditions): {s^ j 
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1. Disorder condition (hot start): In this state we start our program with 
complete disordered state, which means the entropy is maximal. 

2. Order condition (cold start): In this state we start our program with 
complete ordered state, which means the entropy is minimal and the 
temperature is 0. 

2. Generate a new state: 

1. Choose ij th element randomly (Monte Carlo) 

2. Reverse ij th spin direction to create a trial configuration 

3. Calculate the energy of the trial configuration 
4. 
5. Else with random process with P probability accept the change 

(b) Simulation Algorithm [Appendix: Project.f90]: 



If E new < E {si j} then accept the change 




Subroutine Monte Carlo 

Monte Carlo choose of a random particle 
Metropolis to change orientation of particle 



Loop for 
20000000 time 



Spin Structure 




Total energy 
Average magnetization 



Spin_x.dat 



Results.dat 



(c) Metropolis Algorithm [Appendix: Subroutine MonteCarlo]: 



For random point ij: 
s'(iJ) = -s(ij) 




> s(i j) = s'(ij) 



-► loop 
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T=1.0, Average magnetization^ 0.9446 
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T=1.5, Average magnetization^ 0.7316 



T=2.0, Average magnetization^ 0.3589 
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T=4.0, Average magnetization^ 0.1107 
Figure 4. The result structures in 40x40 lattice. 



CONCULSION 



The result structures of simulation (in 40x40 lattice) are shown in Figure 4. As we expect, the 
entropy of system increasing by the temperature of system. 

Figure 5 are the charts of results for 20x20, 40x40, 60x60 and 100x100 in compare of 
analytical graph. As we expect, the value of the average magnetization in high temperature 
are decrease by increasing of system. It's mean that the system becomes nearest to statistical 
limit and the simulation work well as a statistical system. The other point is that the 
qualitative behavior of system and the simulation's results are same. On the other hand, we 
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know from renormalization group method in the critical points we have a self similarity 
property for the system. This figure shows that the critical point is T c = 2.1 j/k B , which is 
near to Curie temperature (T c = 2.269 j/k B ) validly. 
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Figure 5. The result average magnetization simulation results. 
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APPENDIX: Project.f90 

! 

[Computational Physics Project (Fall 2008) 

!By Sayed Mohammad Mahdi Sadrnezhaad - 8711105 

!smmsadr@gmail.com - smmsadr@merc.ac.ir 

! 

!This program writen to calculate the equilibrium states of the Ising Hamiltonian 
!with Monte-Carlo methods. 

I 

module parameters 

!In this module we define the value of our parameters. 

integer N,Total 

real(8) J,Total_energy,kB,T,Average_magnetization 

integer(4) P 

integer s( 100, 100) 
end module parameters 
I 

Program Monte_Carlo_for_Ising_Hamiltoni an 

use parameters 

implicit none 

[integer i 

open (10,file='p.daf,status='old') 
read (10,*) P 

close (10) 

Average_magnetization=0 

Total=0 

open (20,file='Results.dat',status='unknown') 

call parameters_reader 

call Ising_Hamitunian 

close (20) 

open (10,file='p.dat',status='unknown') 
write (10,*) P 

close (10) 
end Program Monte_Carlo_for_Ising_Hamiltonian 



subroutine parameters_reader 
IThis subroutine reads parameters from input file 2. 
use parameters 
implicit none 
integer b,c 
s=0 

open (10,file='parameters.dat',status='old') 
read (10,*) J 
read (10,*) N 
read (10,*) kB 
read (10,*) T 
close (10) 

!open (10,file=First_spin.dat',status='old') 
dob=l,N 

doc=l,N 

!read (10,*) s(b,c) !For when we want read first spin from file. 
s(b,c)=l 
end do 
end do 
!close(10) 
end subroutine 
! 

subroutine Montecarlo (counter) 
use parameters 
implicit none 

integer r_x,r_y, counter !, reject 
real r,Probability,ran 
r_x=int(ran(P)*(N-l)+l) 
r_y=int(ran(P)*(N-l)+l) 
r=ran(P) 
[print* ,N,r_x,r_y 
if (Probability (r_x,r_y)>l) then 

s(r_x,r_y)=- 1 * s(r_x,r_y) 
else 

if (Probability(r_x,r_y)>r) then 

if (Probability(r_x,r_y)<l) s(r_x,r_y)=-l*s(r_x,r_y) 
end if 
end if 

counter=counter+ 1 
[print*, counter 
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end subroutine 

! 



function Probability(x,y) 

use parameters 

implicit none 

integer x,y,xl,yl,x2,y2 

real Probability,delta_E,beta 

beta=-l/(kB*T) 

xl=mod(x+N-2,N)+l 

x2=mod(x,N)+l 

yl=mod(y+N-2,N)+l 

y2=mod(y,N)+l 

delta_E=-J*s(x,y)*(s(xl,yl)+s(xl,y2)+s(x2,yl)+s(x2,y2)) 

Probability =exp(beta* delta_E) 
end function 
j 

subroutine Ising_Hamitunian 
use parameters 
implicit none 

integer counter,b_b_b !,Tc,r_x,r_y,i 
!real r 
b_b_b=100 
Total=0 

Average_magnetization=0 
Total_energy=0 
counter=0 
do while(counter<=20000000) 

call Montecarlo(counter) 
if (counter>b_b_b*1000) then 
b_b_b=b_b_b+l 
call Energy_calculator 
call Magnetization_Calculator [(counter) 
!if (Number_of_Zero>300000) then 
call write_results (counter) 

!Print*,Counter,Average_magnetization,Total_energy 
end if 
end do 
end subroutine 

! 



subroutine Energy _calculator 
use parameters 
implicit none 
integer x,y,xl ,y 1 ,x2,y2 
!real E 

Total_energy=0 
dox=l,N 

doy=l,N 

xl=mod(x+18,N)+l 



x2=mod(x,N)+l 
yl=mod(y+18,N)+l 
y2=mod(y,N)+l 

Total_energy=Total_energy- J* s(x,y) *(s(x 1 ,y l)+s(xl ,y2)+s(x2,y 1 )+s(x2,y2))/4 
end do 
end do 

Total_energy =Total_energy/(N * * 2) 
end subroutine 



subroutine Magnetization_Calculator !(c) 
use parameters 
implicit none 
integer f,g,sum2 !,c 
sum2=0 
dof=l,N 

dog=l,N 

sum2=sum2+s(f,g) 
end do 
end do 

Average_magnetization=sum2+Average_magnetization 
total=total+l 
end subroutine 
j 

subroutine write_results(c) 
!This subroutine writes results in output file, 
use parameters 
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implicit none 
integer i,k,c 

open (30,file='spin_l .dat',status='unknown') 
open (3 l,file='spin_2.dat',status='unknown') 
doi=l,N 

dok=l,N 

if(s(i,k)=l)then 

write (30,4) i,k,s(i,k) 
else 

write (3 1,4) i,k,s(i,k) 
end if 
end do 
end do 
close (30) 
close (31) 

write (20,4) c,Total_energy,Average_magnetization/(total*(N**2)) 
4format(9(2x,gl7.11)) 
end subroutine 

! 



